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ABSTRACT 

The dust properties in the Large and Small Magellanic Clouds are studied using the HERITAGE 
Herschel Key Project photometric data in five bands from 100 to 500 /im. Three simple models of 
dust emission were fit to the observations: a single temperature blackbody modified by a power- 
law emissivity (SMBB), a single temperature blackbody modified by a broken power-law emissivity 
(BEMBB), and two blackbodies with different temperatures, both modified by the same power-law 
emissivity (TTMBB). Using these models we investigate the origin of the submm excess; defined 
as the submillimeter (submm) emission above that expected from SMBB models fit to observations 
< 200/xm. We find that the BEMBB model produces the lowest fit residuals with pixel-averaged 
500 /xm submm excesses of 27% and 43% for the LMC and SMC, respectively. Adopting gas masses 
from previous works, the gas-to-dust ratios calculated from our the fitting results shows that the 
TTMBB fits require significantly more dust than are available even if all the metals present in the 
interstellar medium (ISM) were condensed into dust. This indicates that the submm excess is more 
likely to be due to emissivity variations than a second population of colder dust. We derive integrated 
dust masses of (7.3± 1.7) x 10^ and (8.3dz2.1) x 10^ Mq for the LMC and SMC, respectively. We find 
significant correlations between the submm excess and other dust properties; further work is needed 
to determine the relative contributions of fitting noise and ISM physics to the correlations. 

Subject headings: infrared: galaxies, infrared: ISM, ISM: general, Magellanic Clouds 
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1. INTRODUCTION 

Among nearby galaxies, the Large Magellanic Cloud 
(LMC) and Small Magellanic Cloud (SMC) represent 
unique astrophysical laboratories for interstellar medium 
(ISM) studies. Both Clouds are relatively nearby, the 
LMC at ^50 kpc (Walker 2012) and the SMC at ^60 
kpc (Hilditch et al. 2005), and provide ISM measure- 
ments that are relatively unconfused along the line-of- 
sight as compared to similar observations in the Milky 
Way (MW). The LMC and SMC ultraviolet dust ex- 
tinction properties show strong variations both inter- 
nally and in global averages in comparison to each other 
and the MW (Lequeux et al. 1982; Prevot et al. 1984; 
Clayton & Martin 1985; Fitzpatrick 1985; Gordon et al. 
2003; Maiz Apellaniz & Rubio 2012). The two Clouds 
span an important metallicity range with the LMC at 
^^1/2 Zq (Russell & Dopita 1992) being above and the 
SMC at ^1/5 Zq (Russell & Dopita 1992) being below 
the threshold of 1 /3-1 /4 Zq where the properties of the 
ISM change significantly as traced by the reduction in 
the Polycyclic Aromatic Hydrocarbon (PAH) dust mass 
fractions and (possibly) dust-to-gas ratios (Draine et al. 
2007). The far- infrared (FIR) to submillimeter (submm) 
emission from the Clouds shows more submm emission 
than expected from existing dust grain models, with the 
SMC having a larger amount of this excess emission 
(Israel et al. 2010; Bot et al. 2010a). 

The submm excess was seen first in the MW using 
the COBE/FIRAS (Boggess et al. 1992; Mather et al. 
1993) observations of high-latitude cirrus dust emission 
(Wright et al. 1991; Reach et al. 1995). These works 
found the 100-300 fim observations were well modeled 
with a single temperature blackbody modified with a 
power law emissivity, but that the longer wavelength ob- 
servations (A > 300 /rm) required a second dust compo- 
nent with a temperature of 4-7 K. The spatial correla- 
tion of this second dust component with the hotter main 
dust component along with physical arguments on dust 
heating led Reach et al. (1995) to argue that emissiv- 
ity variations away from a simple power law were more 
likely to explain the observations than a second com- 
ponent of very cold dust. The need for a non-trivial 
FIR to submm dust emissivity shape was quantified by 
Li & Draine (2001) where they modified the emissivity of 
“astronomical” silicate grains to have an emissivity with 
a shallower wavelength dependence at A > 200 ^m than 
at A < 200 fim. More recently, Paradis et al. (2012) an- 
alyzed Herschel Space Observatory (Pilbratt et al. 2010) 
observations of the MW plane and found a significant 
submm excess at 500 /rm that increased from the inner 
to the outer Galaxy. 

Previous work on the submm excess in nearby galax- 
ies by Galliano et al. (2003, 2005) and Galametz et al. 
(2011) used the combination of FIR observations 
(A < 200 /im) from the Infrared Space Observa- 

tory (Kessler et al. 1996) and Spitzer Space Telescope 
(Werner et al. 2004) with submm observations (A ^ 
850 /im) taken using ground-based observatories. These 
works provided strong evidence of a submm excess at 
^850 /im and that this excess is largest in low metal- 
licity galaxies. With the advent of Herschel obser- 
vations, the presence of a submm excess at 500 /im 
has been established in many low metallicity galaxies 


including the Magellanic Clouds (Gordon et al. 2010; 
Meixner et al. 2010; Galliano et al. 2011; Dale et al. 
2012; Kirkpatrick et al. 2013; Remy-Ruyer et al. 2013). 

The dehnition of the submm excess has not been uni- 
formly defined in the literature, complicating the com- 
parisons between different studies. Generally, a model 
is used to define the zero submm excess baseline; this 
model varies from simple modified blackbodies to more 
complex dust grain models. In addition, the uncertain- 
ties assumed on the observations have varied leading to 
the same submm excess level being considered signifi- 
cant by one work and not significant by another. This 
illustrates the need for a uniform definition of reference 
spectral energy distribution (SED) from which to mea- 
sure the submm excess and a common set of assumptions 
on the observational uncertainties. It is also critically 
important to properly include the full observational un- 
certainties, both correlated and uncorrelated, as shown 
by Galliano et al. (2011) and Veneziani et al. (2013). 

For clarity in this paper, we adopt the definition of 
the submm excess as the excess emission seen at submm 
wavelengths above that expected for dust grains with a 
single temperature and a emissivity law. This sim- 
ple model is used to fit an observed SED, with the value 
of /3efi providing a measure of the effective emissivity law. 
The origin of the observed effective emissivity law vari- 
ations may be due to one or a combination of factors 
including intrinsic dust emissivity variations, mixing of 
different dust compositions, and variations in dust tem- 
peratures along the line of sight. 

Laboratory studies of the two main interstellar 
dust analogs have shown that carbonaceous grains 
have /3 ~ 1 — 2 (Mennella et al. 1995; Zubko et al. 
1996; Jager et al. 1998) and silicate grains have 
/3 ~ 2 (Mennella et al. 1998; Boudet et al. 2005; 

Coupeaud et al. 2011) in the EIR and submm wave- 
length range. The value of /3eff for a mixed composi- 
tion dust population is determined by both the actual 
ratio of the two compositions and the spectral shape of 
the heating radiation field. Silicate and carbonaceous 
grains have significantly different ultraviolet/optical ab- 
sorption properties and any change in the radiation field 
spectrum will change the luminosity weighting present 
in the infrared (IR) dust emission SED. Deviations from 
simple emissivity laws and dependence on tem- 

perature are seen in laboratory work on dust analogs, 
with silicate grains having larger such variations than 
carbonaceous grains (Mennella et al. 1998; Boudet et al. 
2005; Coupeaud et al. 2011). Such deviations have al- 
ready been seen in astronomical observations, leading 
Li & Draine (2001) to modify their model of ’’astronom- 
ical” silicates such that it already includes a submm ex- 
cess of 11% at 500 /rm, according to our definition above. 
Similar broken power law dust emissivities have been 
implied by FIR to submm observations of the different 
phases of the MW ISM (Paradis et al. 2009). 

Multiple dust temperatures along the line-of-sight can 
also cause effective emissivity law variations. The sim- 
plest case to consider is two dust populations with the 
second population having a significantly colder temper- 
ature than the first. Fitting the composite SED of 
this dust with a single temperature A“^®® emissivity 
law model will result in a submm excess at the wave- 
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lengths where the second cold dust population con- 
tributes. Such two temperature models have been stud- 
ied by Juvela & Ysard (2012) who find that the /3eff 
can either be higher or lower than the intrinsic j3 de- 
pending the distribution of temperatures. More com- 
plex temperature mixing has been investigated with sim- 
ilar results (Shetty et al. 2009a, b; Juvela & Ysard 2012; 
Ysard et al. 2012). 

The implications for our understanding of dust grain 
properties are quite different depending on the origin of 
the submm excess. If the submm excess is due to very 
cold dust, then the total dust mass would potentially in- 
crease significantly as a large mass of cold dust is needed 
to reproduce the observed emission (e.g., Galliano et al. 
2005). On the other hand, if the submm excess is due to 
dependencies of the effective emissivity law with wave- 
length, then this provides insights into variations in the 
ratio of silicate/carbonaceous grains and/or variations in 
spectral shape of the illuminating radiation field. 

The Magellanic Clouds provide two of the best labora- 
tories to study the submm excess given their proximity 
and lower than MW metallicities. Work on this topic 
in the Magellanic Clouds prior to the Herschel observa- 
tions has used ground-based submm observations (e.g., 
Hot et al. 2010b) or low spatial resolution PLANCK ob- 
servations. In particular, the studies by Israel et al. 
(2010) and Hot et al. (2010a) clearly show a submm ex- 
cess in both Clouds, even though the works were focused 
on the longer wavelength emission of the Clouds. They 
found that the observed submm excess can be explained 
using Draine & Li (2007) models with cold dust grains, 
but not by emission due to spinning grains, which is 
the likely origin of the excess emission they observed 
at millimeter to centimeter wavelengths. Similar re- 
sults for the submm excess in the SMC were found using 
the PLANCK observations (Planck Collaboration et al. 
2011, Verdugo et al. submitted). In apparent con- 
flict with these wide- field and/or global studies of dust 
emission in the Clouds, a spatially resolved study by 
Galametz et al. (2013) found no evidence for a submm 
excess at 870 /rm in N159, a massive star-forming com- 
plex in the LMC. As noted by the authors, however, their 
conclusions apply only to high surface brightness regions 
that can be detected using ground-based submm obser- 
vations. 

The HERschel Inventory of The Agents of Galaxy Evo- 
lution (HERITAGE) in the Magellanic Clouds Herschel 
Key Project has mapped both Clouds providing observa- 
tions at 100, 160, 250, 350, and 500 /xm (Meixner et al. 
2013). The HERITAGE wavelength coverage (100- 
500 /im) and spatial resolution (^10 pc at 500 /xm) is well 
suited to measuring the spatial variations of dust proper- 
ties probed by EIR and submm emission. In particular, 
these observations are ideally suited to investigating the 
nature of the submm excess and how it varies spatially in 
each Cloud. The HERITAGE project test observations 
of a strip in the LMC have been analyzed and a measur- 
able submm excess at 500 /mi was found using both sim- 
ple single temperature blackbodies (Gordon et al. 2010) 
and a more complex dust grain model (Meixner et al. 
2010; Galliano et al. 2011). These studies found that this 
submm excess was anti-correlated with ISM (gas or dust) 
surface density. 

The goal of this paper is to investigate the submm 


excess in both Magellanic Clouds using the full HER- 
ITAGE data using simple dust emission models based on 
one or two modified blackbodies. We choose to use such 
models for this paper since they allow large potential 
variations in the effective emissivity laws, whereas exist- 
ing dust grain models do not incorporate the full range 
of variations indicated by laboratory studies of ISM dust 
analogs. In addition, we are careful to use a robust model 
of the uncertainties in the measurements, including the 
correlations between the different Herschel bands due to 
the absolute flux calibration and the background sub- 
traction. Preliminary versions of the dust surface density 
maps derived in this paper were used to study the corre- 
lation between dust and stellar properties in the Magel- 
lanic Clouds by Skibba et al. (2012). 

2. DATA 

The FIR and submm observations of the Magellanic 
Clouds analyzed in this study were taken as part of the 
HERITAGE Key Project (Meixner et al. 2013) using the 
PACS (Poglitsch et al. 2010) and SPIRE (Griffin et al. 

2010) instruments on the Herschel Space Observatory. 
The observations provided images of the LMC and SMC 
at 100, 160, 250, 350, and 500 /xm that cover the 
entire IR emitting regions of both galaxies (8°x8.5° 
and 5°x5° -I- 4°x3° for the LMC and SMC, respec- 
tively). The observation and data reduction details 
can be found in Meixner et al. (2013). It is useful 
to note that as part of the data reduction, the IRAS 
100 /xm (Schwering & Israel 1989; Schwering 1989) and 
MIPS 160 /xm images (Meixner et al. 2006; Gordon et al. 

2011) for each galaxy were used to correct for the drifting 
baseline of the PACS bolometers. Thus the PACS 100 
and 160 /xm images contain the IRAS 100 and MIPS 160 
information as well as the new PACS observations. 

Additional processing steps were performed for this 
study to create images that had the same spatial res- 
olution and the same foreground/background subtrac- 
tion. First, each image was convolved with a kernel that 
transformed the spatial resolution of the images to the 
lowest resolution of the set of images, set by the SPIRE 
500 /xm point-spread- function (PSF) which has a reso- 
lution of ^40". The Aniano et al. (2011) convolution 
kernels were used for this step as they directly and op- 
timally transform the native PSF to that of the SPIRE 
500 /xm PSF. 

Second, a foreground subtraction was done to remove 
the structured emission due to MW dust (cirrus) emis- 
sion. The detailed structure of the MW dust emis- 
sion in the PACS and SPIRE bands was predicted us- 
ing the integrated MW velocity HI gas maps in the 
direction of the LMC (Staveley-Smith et al. 2003) and 
SMC (Stanimirovic et al. 2000; Muller et al. 2003) and 
the Desert et al. (1990) model for the local interstellar 
radiation field. This model gives the conversion between 
HI column and infrared emission. The conversion coef- 
ficients used were 1.073, 1.848, 1.202, 0.620, and 0.252 
(MJy/sr) (1 X 10^° H I atoms/cm^)“^ for 100, 160, 250, 
350, and 500 /xm, respectively. These conversion coef- 
ficients are higher than those that would be obtained 
with the newer DustEM model (Compiegne et al. 2011) 
for the same radiation field, but are similar to the ob- 
served correlations between the MW velocity integrated 
HI and the diffuse emission measured in the same bands 
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in regions outside of the SMC. This step was particularly 
important for the SMC where structures with similar sur- 
face brightnesses to those in the galaxy were removed by 
this subtraction. 

Finally, residual large scale structure in the back- 
ground was removed using a low order 2D surface poly- 
nomial interpolation that was constrained by regions ex- 
ternal to each galaxy. The baseline subtraction reduction 
step for PACS and SPIRE data used different assump- 
tions for these external regions (Meixner et al. 2013) and, 
thus, this final step ensures that all the images have 
the same background subtraction. This background sub- 
traction is especially important for the LMC where the 
SPIRE observations included emission near the edges 
of the HERITAGE coverage due to the very extended 
nature of the LMC (especially south of the LMC main 
body) and the excellent sensitivity of the SPIRE instru- 
ment. 


3. MODELS 

We use three different models to fit the FIR/submm 
surface brightness measurements. The first model is 
a single temperature blackbody modified by a single 
power law emissivity (SMBB). The second model as- 
sumes the submm excess emission is due to variations 
in the wavelength dependence of the dust emissivity law 
that is parametrized by a broken power law (BEMBB). 
The third model assumes the submm excess emission is 
due to a second, lower temperature population of dust 
grains (TTMBB). All our models assume equilibrium 
heating only and so we restrict our fits to using only 
data >100 /rm. It is reasonable to expect that the emis- 
sion at these wavelengths is dominated by equilibrium 
emission from dust grains. In this analysis, any residual 
100 /im contribution due to emission from transitionally 
heated grains will yield a somewhat higher dust tem- 
perature (and thus a smaller dust column density) than 
would be found with our models. In the great majority 
of sight lines, this contribution is too small to be of con- 
cern, but may introduce a systematic bias in the regions 
near intense star formation. 

In general, the surface brightness of dust with temper- 
ature, Tdi is 


Sx = TxBx[Td) 

(1) 

= NdTTa^QxBx{Td) 

(2) 

= —na^QxBx{Td) 
rud 

(3) 

= 4 Tra^QxBx{Td) 

3^ P 

(4) 

= ^XddQxBx{Td) 
Aap 

(5) 

= K\T,dBx 

(6) 


where t\ is the dust optical depth, Nd is the dust column 
density, a is the grain radius, Q\ is the dust emissivity, 
B\ is the Planck function, is the dust surface mass 
density, rud is the mass of a single dust grain, p is the 
grain density, k\ is the grain absorption cross section per 
unit mass. These equations can be evaluated in stan- 
dards units (e.g. cgs or MKS). We found it convenient 
to express in Mq pc“^, k\ in cm^ g“^, and B\ and 


S\ in MJysr ^ and then Eq. 6 is 

Sx = (2.0891 X lQ-^)nxT.dBx. (7) 

Erom Eq. 6, it is clear that the values of k,\ and are 
completely degenerate. Without further information FIR 
to submm SED observations only constrain t\ = kaE^. 
Breaking this degeneracy is possible in the one envi- 
ronment where we have measurements of the expected 
amount of dust independent from the measured FIR to 
submm dust emission. This environment is the MW dif- 
fuse ISM where ultraviolet and optical gas-phase absorp- 
tion measurements provide a strong constraint on the 
depletions in the ISM (e.g., Jenkins 2009). We use these 
measurements to calibrate k,\ in §5 for the models intro- 
duced below. This calibration ensures that our models 
produce the right E^ in the one place where we know the 
correct value from independent measurements. 


3.1. SMBB: Simple Emissivity Law Model 

The SMBB predicts the surface brightness assuming 
a dust population with single dust temperature modi- 
fied by a simple emissivity law (Hildebrand 1983). The 
adopted emissivity law is 




g 

^eff,160 

160-/3eft 


(8) 


The value of nfg ;^gQ is set by fitting of the diffuse MW 
SED (§5 and Table 2). The full set of fit parameters for 
the SMBB model are 9s = (E^, Teff,d, /3eff)- The values 
for the dust properties are effective values due to com- 
position and temperature mixing along the line-of-sight 
and are not directly comparable to interstellar dust grain 
analogs studied in the laboratory (see §1). 


3.2. BEMBB: Broken Emissivity Law Model 

The BEMBB predicts the surface brightness assuming 
a dust population with a single dust temperature modi- 
fied by a broken emissivity law. The adopted emissivity 
law is 




^eff,160 


E{\) 


(9) 


and 


( 10 ) 


where Xh is the break wavelength and is limited to 
> 175 /im. This emissivity law is similar in form to 
that used by Li & Draine (2001) for astronomical sili- 
cates. The value of igg is set by fitting of the diffuse 
MW SED (§5 and Table 2). 

As we are particularly interested in measuring the 
submm excess, we define the submm excess as the excess 
emission at a particular submm wavelength above or be- 
low that expected for a SMBB model with = /3eff,i- 
Given the BEMBB model definition, the submm excess 
at 500 /rm is 


6500 — 


500 / 


/3eff,2 — /3eff,l 


- 1 . 


( 11 ) 


Using 6500 as one of the fit parameters (instead of /3eff,2), 
the fit parameters for the BEMBB model are 6 *be = 
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(Ed,Teff,d,/3eff,i, Ah,e5oo). Note that the value of 6500 
can be negative and this would indicate a submm deficit. 
The values for the dust properties are effective values due 
to composition and temperature mixing along the line- 
of-sight and are not directly comparable to interstellar 
dust grain analogs studied in the laboratory (see §1). 


3.3. TTMBB: Two-Temperature Model 

The TTMBB predicts the surface brightness assum- 
ing two dust populations with distinctly different dust 
temperatures modified by a single, non-broken emissiv- 
ity law. The surface brightness is then 

Sx = Ka [EdiHA(7eff,dl) + Ed2SA(Teff,d2)] (12) 


where 


KA 


^eff,160 ' 
160-Aft 


(13) 


the subscripts d\ and d2 refer to the two dust compo- 
nents, and Teff^di > Teff_d2- The value of is set by 

fitting of the diffuse MW SED (§5 and Table 2). 

Eor this model the submm excess at 500 /im is 


6500 — 


Erf2l?50o(Teff,d2) 

^dlB^OoiTeS^dl) 


(14) 


Again, we use 6500 as a fit parameter and the full set 
of fit parameters for the TTMBB model are 0 tt = 
(Edi, Teff,d2, l3eS, 6500)- Note that the value of 6500 

for the TTMBB model cannot be negative unlike the case 
for the BEMBB model. The values for the dust proper- 
ties are effective values due to composition and temper- 
ature mixing along the line-of-sight and are not directly 
comparable to interstellar dust grain analogs studied in 
the laboratory (see §1). 


3.4. Restricted Models 

It is often assumed in modified blackbody fitting that 
only /3eff values between 1 and 2 are valid. This is 
based on arguments that laboratory measurements of 
dust analogs only give /3 values between these limits. 
More precisely, laboratory measurements of carbona- 
ceous and silicate dust analogs give j3 values between 
0.8 and 2.5 for the Herschel wavelength range (e.g., 
Jager et al. 1998; Coupeaud et al. 2011). It is clear 
that luminosity weighted mixing of dust analogs with 
P values between 0.8 and 2.5 will always result in /3eff 
values in the same range. Yet this is not necessarily 
the case for temperature mixing along the line-of-sight 
(Juvela & Ysard 2012). Combining the effects of compo- 
sition and temperature mixing using full radiative trans- 
fer models, Ysard et al. (2012) give evidence that find 
that an /3eff (/3coior in their terminology) between 0.8 and 
2.5 is reasonable for a range of realistic cases. Thus, we 
include versions of the SMBB, BEMBB, and TTMBB 
models that have /3eff values restricted to be between 0.8 
and 2.5. But we caution that it is more statistically cor- 
rect to include /3eff values outside this range as measure- 
ment noise can create SEDs that require non-physical 
/3eff values to provide statistically robust fits. 


3.5. Band Integration 

Our models produce SEDs that are well sampled in 
wavelength, but our observations have a very coarse 


Table 1 
Grid Parameters 


Parameter 

Range 

Spacing 

SMBB 

log(Sd) [M© pc P 

-4 to 1 

0.1 

Tefl.d [K] 

5 to 75 

1 

/3eff 

-1 to 4 

0.25 

BEMBB 

log(Sd) [M© pc P 

-4 to 1 

0.1 

Tefl,d [K] 

5 to 75 

1 

PeS,l 

-1 to 4 

0.25 

\b [pm] 

175 to 375 

25 

6500 

-1 to 2 

0.25 

TTMBB 

log(Srfi) [M© pc P 

-4 to 1 

0.1 

Teff.dl [K] 

5 to 75 

2 

Teff.dl [K] 

4 to 75 

2 

deff 

-1 to 4 

0.25 

6500 

0 to 2 

0.25 


wavelength sampling as they are taken through filters 
with broad response functions. It is important to cor- 
rectly model the effects of these broad response func- 
tions on the models to give accurate fits to the observa- 
tions. For this paper, we start with the model predic- 
tions of the surface brightnesses at a wavelength resolu- 
tion that well resolves the BAGS and SPIRE bandpasses 
(Muller et al. 2011b; Griffin et al. 2013). Then, the band 
surface brightnesses were determined by integrating over 
their respective band response functions using 


/ S„RE{i^)du 
~ K^olv)-^RE(y)dv 


(15) 


where Re{^) is the response function appropriate for 
extended sources given in fractional transmitted energy. 
The Eo = c/Xo values are given by Ao = 100, 160, 250, 
350, and 500 ^m for the bands with the same names. 
Eq. 15 mathematically models the data that is produced 
by the PACS and SPIRE instruments and data reduc- 
tion pipelines. The integration is done in energy units 
(e.g., MJy sr“^) as both instruments use bolometers that 
measure energy (not photons). The denominator of this 
equation normalizes Re{X) and accounts for the PACS 
and SPIRE calibration convention where the calibration 
is given at specific wavelengths (Aq) and for a S{u) = 
reference spectrum. 


4. FITTING TECHNIQUE 

We computed the models on discrete grids with spac- 
ings fine enough to resolve the final ID likelihoods for 
each parameter. The grids were computed over a large 
range in each parameter to ensure that the likelihood 
function was well sampled. The ranges and spacings for 
both models are given in Table 1. We use a logarith- 
mic spacing for E^ to provide a computationally efficient 
sampling of the full dynamic range of this parameter. 
The minimum and maximum ranges of the parameters 
were set iteratively, expanding the fit parameter ranges 
until the ID likelihood function for the vast majority of 
the pixels in the galaxies were well sampled. 

We fit each pixel that was detected at 3cr above the 
background in all five bands. The probability that a 
particular model fits the data was computed assuming a 
multi-variate Normal/Gaussian distribution (Gut 2009) 
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using 


where 


and 


1 


1 


P(S°^^\0) = ^exp(--x^(0)), (16) 


= (27r)’"det |C| (17) 


/(0) = - 5“°‘^(0)]. (18) 


^obs jg observed SED for a single pixel in the n = 5 

bands, S™°‘^ is the SED for a particular model and pa- 
rameter set, 9, and C is the covariance matrix. The ^ 
notation denotes the transpose of the vector. The covari- 
ance matrix is often given as the S symbol, but we have 
chosen to use C to avoid confusion with the dust surface 
density or standard summation symbol. 

The explicit use of a covariance matrix in the fitting al- 
lows us to directly account for correlations between bands 
in the data. This is a different approach than has been 
recently taken by other authors. One technique for in- 
vestigating the effects of correlated noise on model fit 
parameters is to perform many Monte Carlo trials of 
the observations where they are perturbed by the ran- 
dom and correlated noise and fit with the model (e.g., 
Galliano et al. 2011). A second technique is to include 
parameters in a hierarchical Bayesian model for the cor- 
relations in the absolute flux calibration between bands 
and then marginalize (integrate) over them to deter- 
mine their final fit probabilities (e.g., Kelly et al. 2012). 
While not often done, it is critical to account for the 
correlated noise in observations as neglecting such noise 
terms can significantly bias the resulting fit parameters 
(Veneziani et al. 2013). By including the covariance di- 
rectly into the likelihood function we do not need to per- 
form many Monte Carlo trials for every pixel or use a hi- 
erarchical Bayesian model to account for this noise term. 
In other words, we can include the correlations directly 
in the individual fits efficiently without having to appeal 
to the ensemble behavior. 


2011a; Balog et al. 2013). Similar to SPIRE, for ex- 
tended sources we add an additional 5% correlated un- 
certainty to account for uncertainties in the total beam 
area resulting in a 10% correlated uncertainty between 
bands. Einally, we assume the PACS and SPIRE cali- 
brations are independent given that PACS is calibrated 
using stars and SPIRE using Neptune. Given this infor- 
mation the elements of Ccai are 


(Ccal),, = Sr^{9)S';;^°^{9) [(Aeor)y + (Aancor)*,] (20) 
where 



■ 0.l2 

0.l2 

0 

0 

0 


0.l2 

0.l2 

0 

0 

0 

A — 

■^cor — 

0 

0 

0.082 

0.082 

0.08 


0 

0 

0.082 

0.082 

0.08 


0 

0 

0.082 

0.082 

0.08 


and 


A 


uncor 


■ 0 . 02 ^ 0 0 0 0 

0 0 . 02 ^ 0 0 0 

0 0 0.0152 0 0 

0 0 0 0.0152 0 

0 0 0 0 0.0152 


• ( 22 ) 


The background covariance matrix, Cbkg is calculated 
empirically from a large set of pixels visually identified 
as lying outside of the emitting region of each galaxy. 
The background pixels are in the full images and were 
processed as described in §2. The terms of the covariance 
matrix are calculated using 

. Ef (s? - (&» (si - (Si)) 

where N is the number of background pixels, Sf/S^ is 
the ith/jth band of the /cth pixel, and (Si)/{Sj) is the 
average background in the fth/7'th band. Eor the LMC, 
N = 52113 and 


4.1. LMC and SMC Covariance Matrices 
Eor this work, the covariance matrix is defined as 

C = Ccal + Cbkg (19) 

where Ccai is the absolute surface brightness covariance 
matrix and Cbkg is the background covariance matrix. 
The units of these covariance matrices are (MJy/sr)2. 

The Ccai is given by the details of the PACS and 
SPIRE absolute flux calibrations. The SPIRE instru- 
ment has been calibrated using a model of Neptune 
with an absolute uncertainty correlated between bands 
for point sources of 4% and a repeatability that is un- 
correlated between bands of 1.5% (Griffin et al. 2013; 
Bendo et al. 2013). For extended sources, it is rec- 
ommended to add an additional 4% to account for 
the correlated uncertainty in the total beam area re- 
sulting in an 8% correlated uncertainty between bands 
(Herschel Space Observatory 2011). The PACS instru- 
ment has been calibrated using models of stars and 
asteroids with an absolute uncertainty correlated be- 
tween bands for point sources of 5% and a repeata- 
bility uncorrelated between bands of 2% (Muller et al. 


Cbkg (CMC) = 


and for the SMC, 


Cbkg(^MC) = 


r4.23 

0.78 

0.65 

0.33 

0.14 -| 

0.78 

2.37 

0.85 

0.43 

0.18 

0.65 

0.85 

0.91 

0.47 

0.20 

0.33 

0.43 

0.47 

0.25 

0.11 

.0.14 

0.18 

0.20 

0.11 

0.057. 


N = 4012 and 


r 2.64 

0.56 

0.30 

0.14 

0.064 

0.56 

1.18 

0.46 

0.23 

0.094 

0.30 

0.46 

0.36 

0.20 

0.089 

0.14 

0.23 

0.20 

0.12 

0.054 

. 0.064 

0.094 

0.089 

0.054 

0.030 


(24) 


(25) 

These empirical covariance matrices illustrate that back- 
ground is highly correlated with the correlation increas- 
ing in strength towards longer wavelengths. This is illus- 
trated by the correlation matrix (terms are Cij /[uiUj]) 
for the SMC: 


corrbkg(<S'MC') 


r 1.00 

0.32 

0.31 

0.25 

0.23-1 

0.31 

1.00 

0.70 

0.61 

0.49 

0.30 

0.70 

1.00 

0.94 

0.85 

0.25 

0.61 

0.94 

1.00 

0.91 

. 0.23 

0.49 

0.85 

0.91 

1.00. 
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The EMC correlation matrix is very similar and so is not 
shown. The positive and non-zero correlation terms is 
a signature that the correlated noise in the background 
is due to real astronomical signals. In this case, it is 
traceable to the residual foreground MW cirrus emission 
and the integrated emission from background galaxies. 
The higher covariance values for the EMC is a reflection 
of the increased difficulty of background subtraction for 
this galaxy. 

4 . 2 . Example Fitting Results 

The fitting technique we use fully computes the uD 
likelihood function that a particular model fits the SED 
of a pixel where n is the number of ht parameters. One 
way to visualize the results is to create ID likelihood 
functions for each ht parameter by marginalizing (inte- 
grating) over all the other parameters. This is shown in 
Fig. 1 for the BEMBB model for a single pixel in the 
SMC for three different assumptions; assuming uncor- 
related uncertainties, including the full covariance, and 
including the full covariance while restricting /3eff,i to 
vary between 0.8 and 2.5. The results for pixels in the 
EMC are similar. With the same overall uncertainties, 
we obtain a much narrower function with a stronger like- 
lihood by including the known covariance between the 
bands (§4.1) than by assuming that there is no correla- 
tion between bands. In this case, including the known 
covariance between bands results in better constraints on 
the ht parameters as the allowed model space is reduced. 
The impact of a limited /3eff,i range is shown in this hg- 
ure where, not surprisingly, it makes for a narrower ID 
likelihood function than allowing /3eff,i to vary to fully 
sample the /3eff,i ID likelihood function. Note that this 
limitation simply crops the /3eff,i ID likelihood function, 
but changes the shape of the other ID likelihood func- 
tions signihcantly. 

4 . 3 . Sensitivity Tests 

The goal of the sensitivity tests is to determine if there 
are systematic shifts in recovered parameters and if the 
uncertainty on the recovered parameters matches that 
measured from the widths of the ID likelihood functions. 
We simulated observations by picking a model SED and 
adding noise using the Cholesky factorization of the co- 
variance matrix appropriate as if the model was observed 
like the SMC was observed in HERITAGE. The results 
using the EMC noise model give very similar results. We 
repeated the simulation for each model SED 20 times 
to provide a good sampling of the recovered fit parame- 
ter uncertainties and systematic offset from the input fit 
parameters. 

As we are testing the ability of this fitting technique 
to recover parameters by fitting simulated observations, 
this requires a way to measure the recovery of the input 
model parameters. The main output of the fitting is the 
nD likelihood function, but it is often useful to distill 
these results to “best fit” or summary values. We use 
three different ways to define the “best fit” values. The 
first is the most traditional definition of the “best fit” and 
corresponds to the maximum likelihood (‘max’). This is 
also called the “traditional method in some papers 
(e.g., Kelly et al. 2012; Juvela et al. 2013). The ‘max’ 
value is most useful when plotting the best htting model 


with observations or investigating the fitting residuals. 
The second is the expectation value (‘exp’) which is the 
likelihood weighted average of the parameter and is a re- 
flection of the full likelihood function. This ‘exp’ value 
reflects the best “average” value as it reflects the full like- 
lihood function (not just the peak like the ‘max’ value). 
We End the ‘exp’ particularly useful for making images 
of the fit parameters. The third way to reflect the best 
fit is take a realization of the full nD likelihood function 
itself (‘realize’). This involves randomly sampling the 
likelihood function and reflects the full likelihood func- 
tion’s shape in a statistical sense. The ‘realize’ method 
is most useful when studying the ensemble behavior of 
the fit parameters for many pixels. 

The results for runs with 2000 randomly picked BE- 
MBB models are shown in Fig. 2. All three different 
methods of determining the “best fit” parameters give 
similar results with similar trends with each parameter. 
The ‘exp’ gives the lowest systematic error in the recov- 
ery, but the ‘max’ gives the lowest scatter. The ‘realize’ 
method provides a nominally worse recovery than both 
the other methods, but is a fuller picture of true sen- 
sitivity of the fitting. Overall, which “best fit” method 
used depends on the particular question being asked. We 
illustrate this later in this paper and in the companion 
paper on the gas-to-dust ratio (Roman-Duval et. ah, this 
issue). 

Of particular interest for this paper is the fact that the 
recovery of the submm excess, esoo, is good to around 
10%, on average, for the ‘realize’ method and around 1% 
for the ‘exp’ method. For the companion paper (Roman- 
Duval et al, this issue), the fit parameter of main interest 
is Ed and the recovery is good, on average and in log(Ed) 
units, to 0.05 for the ‘realize’ method and 0.001 for the 
‘exp’ method. This excellent recovery of log (Ed) holds 
even in the presence of significant scatter in Teff.d and 
may be due to other parameters in the fitting varying 
to compensate. Note that for the ‘exp’ method we com- 
puted the expectation value of log(Ed) as we found that 
the sensitivity tests showed significantly less systematic 
bias than if we computed the expectation value of Ed. 
We confirmed that the widths of the ID likelihood func- 
tions matches the noise in the recovery of the input model 
parameters. 

4 . 4 . Number of Parameters and Data Points 

The number of parameters in our models is three, five, 
and five for the SMBB, BEMBB, and TTMBB models, 
respectively. In this paper, these models are fit to FIR- 
submm SEDs that are composed of five data points. At 
first glance, this violates the rule that fitting requires at 
least one data point more than the number of fitting pa- 
rameters to provide a unique solution. This is correct, 
if the fitting is done with a model that can fit any dis- 
tribution of data points. This is clearly not the case for 
our models as they are all constrained to have a spec- 
tral shape of one or two modified blackbodies. In other 
words, they cannot fit arbitrary spectral shapes but are 
constrained by our knowledge of the physics of dust grain 
emission. Effectively, we are using more than just five 
data points in our fits as we combine the data points 
with a larger body of observations that informs our un- 
derstanding of dust physics and, therefore, the appro- 
priate models to use. Finally, our use of full likelihood 
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Figure 1. The ID likelihood functions for a single pixel in the SMC using the BEMBB model are plotted for fitting while assuming 
uncorrelated uncertainties, including the full covariance, and including the full covariance while restricting the allowed /3efl 1 values to be 
between 0.8 and 2.5. Note that /?efF 2 is completely determined by the value of /3eff i and esoo and we present the fSgg 2 ID likelihood 
function for completeness. 


functions explicitly accounts for the impact of the num- 
ber of parameters on how well we can determine each fit 
parameter. Using full likelihood functions has the addi- 
tional benefit of measuring how well each parameter is 
constrained by the data explicitly. Some parameters are 
better constrained than others as shown in Fig. 1. For 
example, and Teff.d are better constrained as the over- 
all level and spectral shape are well constrained by the 
observations, but the detailed spectral shape is less well 
constrained and this impacts / 3 eff,i, Ah, and 6500 strongly. 

5. MODEL CALIBRATION 

It is important to calibrate dust models to repro- 
duce observations where there are independent measure- 
ments of the same quantities using the same fitting tech- 
nique. This is regularly done when setting up full dust 
grain models (e.g., Li & Draine 2001; Zubko et al. 2004; 
Compiegne et al. 2011). One key calibration source is 
the FIR-submm SED of the MW diffuse ISM. This is 
a unique environment as it is the one place where the 
amount of dust has been measured using ultraviolet and 
optical gas-phase absorption lines and knowledge of the 


total amount of atoms expected in the ISM (e.g., Jenkins 
2009). Thus, fitting the FIR-submm MW diffuse SED 
results in a calibration of the dust emissivity k \ as the 
degeneracy between this quantity and is removed. 

In full dust grain models, the calibration of k \ is usu- 
ally set such that the luminosity weighted average re- 
sponse of the different dust grain components reproduces 
the MW diffuse SED when the dust is illuminated by the 
average MW radiation field. In a similar manner, the 
«ieff,i 60 for the models used in this paper is set such that 
htting the MW diffuse SED produces the observed gas- 
to-dust ratio. By determining Keff,i 60 using the measure- 
ments of the diffuse MW emission for each of our models, 
we ensure that our models derive the correct dust surface 
density in the one physical environment where we have 
independent constraints on the dust mass. It is critical to 
note that this calibration does not impose a gas-to-dust 
ratio calibration on our model, just a calibration that we 
derive the correct mass of dust in the MW diffuse ISM. 

This calibration does mean that we are assuming that 
the dust properties in the Magellanic Clouds are the 
same as those in the diffuse MW. This assumption is 
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Figure 2. The results for sensitivity tests of the BEMBB model for 2000 models randomly selected from the full model grid are shown. 
The results are plotted as averages and standard deviations of the recovered minus input parameters in 10 bins over the parameter range. 
The three different methods of determining the accuracy of the recovered parameters are ’max’ = maximum likelihood, ’exp’ = expectation 
value, and ’realize’ = one realization based on the ID likelihood functions for each parameter. 
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reasonable given the evidence from ultraviolet extinc- 
tion measurements in all three galaxies. The SMC does 
show UV extinction curves most different from the av- 
erage in the MW, but it also has curves that are very 
similar to the MW average (Gordon & Clayton 1998; 
Maiz Apellaniz & Rubio 2012). The LMC shows extinc- 
tion curves that are similar or equivalent to the MW 
average (Misselt et al. 1999; Gordon et al. 2003). While 
many of the MW lines-of-sight show extinction curves 
similar to the MW average by definition (Valencic et al. 
2004), there is one line-of-sight that shows a UV extinc- 
tion curve indistinguishable from the most different SMC 
extinction curves (Valencic et al. 2003). It is not clear 
if the globally average UV dust extinction is different 
between the three galaxies, mainly due to small sam- 
ples sizes of such measurements in the Magellanic Clouds 
(Gordon et al. 2003). One piece of evidence that far-IR 
emissivity of dust grains is similar between the MW and 
SMC is the similarity of their Keff.ieo values as derived 
using dust grain model fitting (see §5.3). While it is rea- 
sonable to assume the dust is similar in all three galaxies, 
it is an assumption and the dust surface densities will 
vary inversely in direct proportion to any changes in the 
adopted /Ceff.ieo calibration. 

Evidence for different than the MW dust in the 
LMC was found in work by Meixner et al. (2010) and 
Galliano et al. (2011) using the HERITAGE test obser- 
vations of a strip in the LMC. These works used two mod- 
els of dust, one composed of silicates, graphite, and PAH 
grains that describes average MW dust ( “standard” ) and 
a second with amorphous carbon instead of graphite 
(“AC”). The analysis found that the gas-to-dust ratio 
for the “standard” model was lower than reasonable for 
the LMC metallicity, while the “AC” model produced a 
reasonable ratio. We discuss the issue of gas-to-dust ra- 
tios for the LMC and SMC using the fitting results for 
the models used in this paper and calibrated using the 
MW diffuse SED in §6.3. In addition, we have estimated 
the systematic error on Keffpeo due to assuming that the 
dust is like that in the MW in §5.3. 

Direct measurements of ISM depletions in the Mag- 
ellanic Clouds would allow us to directly calibrate our 
models in these galaxies. This would remove the as- 
sumption that the dust grain compositions in the Mag- 
ellanic Clouds are the same as those in the Milky 
Way. Currently, there exists only a limited number 
of sightlines and atoms with measured depletions in 
the Magellanic Clouds (Roth & Blades 1997; Welty et al. 
1997, 2001; Sofia et al. 2006; Peimbert & Peimbert 2010; 
Welty & Crowther 2010). Extending these studies in 
terms of atomic species and galactic environments should 
be a priority for the astronomical community, since they 
are critical for interpreting the wealth of FIR to submm 
ISM observations obtained by recent space missions. 

5.1. Milky Way Diffuse SED 

For the diffuse MW emission, we use the 
Compiegne et al. (2011) measurement where emis- 
sion was measured by correlating the IR versus HI 
emission maps in atomic gas dominated regions of the 
MW. The IR measurements we use are mainly the 
COBE/FIRAS spectrophotometry from 127 to 1200 /rm 
supplemented by the DIRBE 100 fj,m photometry. 
As we want to calibrate our models using the same 



X [/^m] 


Figure 3. The observed MW diffuse SED from COBE FIRAS and 
DIRBE is plotted along with the best fits for the models used in this 
paper. The best fit is defined using the ‘max’ method discussed in 
§4.3. The ’PACS/SPIRE phot.’ points (purple squares) are those 
used to constrain the fits of the models and were derived from the 
COBE FIRAS and DIRBE measurements. 

bands as used for the HERITAGE observations, we 
integrated this diffuse MW SED using the method 
described in 3.5 for all the bands except the PACS 
100 fj,m band. For this band, we adopted the DIRBE 
100 /im measurement as the bandpasses are similar. 
The resulting MW diffuse SED is 0.71, 1.53, 1.08, 0.56, 
and 0.25 MJy sr“^ (10^° H atom)“^ for the 100, 160, 
250, 350, and 500 ^m and is plotted in Fig. 3. These 
values differ from those given for the same bands by 
Compiegne et al. (2011) mostly as we have not included 
the 0.77 correction for ionized gas. In addition, there 
are minor differences in the response curves used. We 
do not include the 0.77 correction for ionized gas as the 
depletion measurements do not include any ionized gas 
correction. For the uncertainties, we have assumed a 
5% correlated and a 2.5% uncorrelated terms (see §3.5) 
given the high quality of the COBE FIRAS and DIRBE 
calibrations. 

5.2. Milky Way Diffuse Gas-to-Dust Ratio 

As the MW diffuse SED is measured as a correlation 
between dust and gas emission, the constraint we need 
is the MW diffuse gas-to-dust ratio. We use the work 
of Jenkins (2009) to determine the appropriate gas-to- 
dust ratio since this work provides an excellent compila- 
tion and summary of MW depletions. The observed H 
columns of our adopted FIR-submm MW diffuse SED are 
log[fV(iJ)] < 20.7. The average depletion of all the sight- 
lines with these column densities tabulated by Jenkins 
(2009) is F* = 0.36. F* is the depletion factor and mea- 
sures the overall depletions in a sightline. Using the de- 
pletion fits of Jenkins (2009) with F* = 0.36, the diffuse 
MW gas-to-dust ratio is computed to be 150. 

5.3. Calibrating Keffneo 

We calibrate the value of Keff.ieo in each of our models 
so that they reproduce the MW diffuse observed gas-to- 
dust ratio of 150. For our work, we have chosen 160 /rm 
to set our normalization of Keff,A as shorter wavelengths 
have a weaker dependence on temperature based on lab- 
oratory investigations of dust analogs (Coupeaud et al. 
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Table 2 

MW Diffuse Fit Results 


Model 

^efF,160^ 
[cm^ g~^] 

Other Parameters 

Expectation Values 

SMBB 

BEMBB 

TTMBB 

TTMBB 

9.6 ±0.4 ±2.5 

11.6 ±1.5 ±2.5 
517 ±214 ±2.5 

9.6 ±0.4 ±2.5 

('Cff.d, ^efl) 

(T’eff ,d , deff ,1 Ab , 6500 ) 
(T'eff.dl , Teff,d2, /3eff I 6500) 

adopted 

(17.2 ±0.4 K, 1.96 ±0.10) 

(16.8 ± 0.6 K, 2.27 ± 0.15, 294 ± 29 Atm, 0.48 ± 0.11) 
(15.0 ±0.7 K, 6.0 ±0.8 K, 2.9 ± 0.1, 0.91 ± 0.25) 


The results are given as value ± fitting uncertainty ± systematic uncertainty 


2011). The Keff,i 60 values required for each model based 
on the ‘exp’ method of determining the best fits (see 
§4.3) are given in Table 2. The second uncertainty on 
tteff,i 60 is an estimate of the systematic uncertainty (see 
next paragraph). The fit parameters for each model are 
also given in this table, along with Icr uncertainties. The 
larger relative uncertainties on Keffpeo for the BEMBB 
model as compared to the SMBB can be directly traced 
to the larger number of BEMBB ht parameters. The 
‘max’ best fit models are plotted in Fig. 3. 

The /teff.160 values for the SMBB and BEMBB mod- 
els agree favorably with other determinations while the 
value for the TTMBB model does not. For example, 
if “astronomical” silicate grains with a = 0.1 /im and 
p = 3 g cm“^ are used, then Keff.ieo = 13.75 cm^ g“^. 
Such grain properties are often assumed for simple mod- 
ified blackbody hts as this is the average size for a 
Mathis et al. (1977) grain size distribution (Hildebrand 
1983). The widely used Weingartner & Draine (2001) 
full dust grain model for R(V) = 3.1 has a Keff.ieo = 
9.97 cm^ g“^. The updated version of this model has a 
fi^eff,i 60 = 12.5 cm^ g“^ (Draine & Li 2007; Draine et al. 
2014). The Keff.ieo values for the Zubko et al. (2004) 
models that include graphite and amorphous carbon 
range from 10.75 to 15.0 cm^ g“^. Finally, the 
Weingartner & Draine (2001) model for the SMC Bar 
extinction curve with no 2175 A extinction feature has 
fi^eff .160 = 13.1 cm^ g“^. Using the range of these model 
fi^eff,i 60 values we estimate that there is a ±2.5 cm^ g“^ 
additional uncertainty on Keff,i 60 due to systematic un- 
certainties in our knowledge of dust grains. 

The TTMBB model with Keff.ieo = 517 ± 214 cm^ g“^ 
requires a dust grain that is very efficient at emission, 
yet this level of efficiency is much higher than any as- 
tronomically reasonable dust grain. A much simpler ex- 
planation is that the dust in the MW diffuse ISM is not 
well modeled by a TTMBB model that includes a very 
cold (Teff^d ~ 6 K) dust grain population. This is the 
same conclusion given by the Reach et al. (1995) analy- 
sis of the FIRAS data. There still may be regions in the 
ISM of the MW or other galaxies that are well described 
by the TTMBB model. To allow for such regions, we 
adopt the Keff,i 60 of the SMBB model as the value for 
the TTMBB model. 

The variations in the Keffpeo values in the literature 
and between the different models used in this paper 
clearly indicate that /teffpeo is sensitive to the model 
assumptions. Thus, it is important to calibrate each 
model explicitly with the diffuse MW SED and a de- 
pletion measured gas-to-dust ratio. This is a standard 
calibration method for dust grain models (Draine & Li 
2007; Compiegne et al. 2011) and we advocate that such 


calibrations be done for all dust emission models (Bianchi 
2013). Such model calibrations will allow for meaningful 
comparisons between the results from different models. 

6. RESULTS 
6.1. Fitting Residuals 

One obvious question is: Which model, SMBB, BE- 
MBB, or TTMBB, fits the observations best? The an- 
swer to this question will give an indication of the ori- 
gin of the submm excess. The most straightforward 
method to test how well a model fits the data is to 
examine the residuals of the data to the fits. The 
value computed using eq. 18 gives such a quantitative 
measure of the residuals. For the SMC, the pixel aver- 
aged value is 3.47 for the SMBB model, 0.88 for the 
BEMBB model, and 1.83 for the TTMBB. The models 
with 0.8 < /3eff < 2.5 have higher average values than 
the unconstrained versions. For example, the PeB con- 
strained version of the BEMBB model for the SMC has 
an average x^ value of 1.32. The LMC average x^ values 
behave similarly. 

More evidence that the BEMBB fits the data best (out 
of the three models) can be found by examining the be- 
havior of the fit residuals versus surface brightness. Fig. 4 
shows the fit residuals for the SPIRE 250 pm band for 
all three models used in this paper for both Magellanic 
Clouds. The trends for other bands are similar, espe- 
cially in the relative behavior of the fit residuals be- 
tween the models. This figure clearly shows that the sim- 
plest model (SMBB) has residuals larger than expected 
given the known uncertainties. This holds for PeS un- 
constrained and constrained to be between 0.8 and 2.5. 
In addition, the residuals for the SMC have a systematic 
trend with more negative residuals at intermediate sur- 
face brightnesses. Such a trend is not consistent with the 
uncertainties in the absolute flux calibration or the back- 
ground subtraction. Of all models, the BEMBB model 
without any constraint on /?eff fits the data best. Over- 
all, the BEMBB model shows the smallest residuals with 
no obvious trend with surface brightness unlike the other 
models. The BEMBB model consistently shows smaller 
residuals in all the bands, not just the SPIRE 250 pm 
band. The other models have higher overall residuals 
and show systematic offsets and/or trends with surface 
brightness. The BEMBB and TTMBB models have the 
same number of fit parameters, yet the behavior of their 
residuals are different. This illustrates that it is not only 
the number of ht parameters that is critical for the htting 
accuracy, but the allowed spectral shapes. 

Overall, the BEMBB spectral shapes ht the data better 
than the TTMBB and SMBB spectral shapes. This is 
evidence that the submm excess is more likely to be due 
to emissivity variations than a second population of cold 
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Figure 4. The fractional residuals for the SMC (top) and LMC (bottom) of the fits for the SPIRE 250 fim band are shown for all the 
models. Each model has been plotted shifted by multiples of 0.5 on the x-axis. The false color gives the log density of points and each 
point represents the residual for the ‘max’ estimator for a single pixel. The ‘max’ estimator was used to give each model the best chance 
to have the lowest residuals. The plots at other wavelengths show similar behaviors with the BEMBB model having the lowest residuals. 


dust. 

6.2. Total Dust Masses 

The total dust masses are of interest for studies of 
the lifecycle of dust in the LMC and SMC (Boyer et al. 
2012; Matsuura et al. 2013; Zhukovska & Henning 2013). 
In addition, they can be used along with the total gas 
masses as a way to tell if a model produces realistic 
amounts of dust (see §6.3). 

We give the dust masses for the different models in Ta- 


ble 3 integrated over the >3cr pixels. The restricted /3eff 
version of the models produces results that are very sim- 
ilar and are not given in the table. The dust mass values 
are given as total ± statistical uncertainty ± uncertainty 
due to the Keff.ieo uncertainty. To convert from dust 
surface density to dust mass we use distances of 60 kpc 
(Hilditch et al. 2005) and 50 kpc (Walker 2012) for the 
SMC and LMC, respectively. The total dust masses are 
computed from the ‘realize’ method to produce dust sur- 
face density maps that provide a full accounting of the 
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Table 3 

Integrated Dust Masses and Gas-to-Dust Ratios 
Integrated over >3 it Pixels 


Model 

Md [Mq] 

Gas/Dust® 

LMG 

SMBB 

(8.1 ±0.07 ±2.1) X 10® 

340 ± 90 

BEMBBt' 

(6.7 ± 0.03 ± 1.7) x 10^ 

400 ± 100 

TTMBB 

(1.2 ±0.01 ±0.3) X 10'^ 

22 ±6 

expected: scaling MW gas-to-dust ratios 
expected: MW depletions and LMC abundances 
expected: all metals in dust 

200-500 

150-360 

>105 

SMG 

SMBB 

(8.1 ±0.1 ±2.1) X lO'* 

1440 ± 380 

BEMBB*’ 

(6.7 ±0.1 ± 1.7) X 10^ 

1740 ± 440 

TTMBB 

(5.1 ±0.3 ±1.3) X 10** 

230 ± 60 

expected: scaling MW gas-to-dust ratios 
expected: MW depletions and SMC abundances 
expected: all metals in dust 

500-1250 

540-1300 

>300 


The integrated gas masses in Mq for the same areas and 
with the same background removal in the LMC/SMC are 2.5 X 
10®/1.0 X 10® for HI and 2.1 X 10’^/1.6 X 10^ for H 2 (Leroy et al. 
2007a; Hughes et al. 2010). 

Model favored from the analysis in this paper (see §6.1 and §6.3) 

likelihood functions for all pixels. Ten different maps 
were made for each galaxy using the ‘realize’ method that 
samples the likelihood function once for each pixel. This 
provides a robust measurement of the impact of the ht- 
ting noise of each pixel in the integrated dust mass mea- 
surement. The average and statistical uncertainty of the 
integrated dust mass were computed from the ten maps. 
The large number of pixels in each galaxy results in the 
total dust mass changing only slightly between different 
realizations and this is the origin of the small statistical 
uncertainty. These dust masses are integrated only over 
the areas that were detected at 3a above the background 
in ah Hve Herschel bands measured by HERITAGE. Pix- 
els above >3cr contribute 0.79, 0.73, 0.62, 0.61, and 0.61 
of the SMG global fluxes of 15.7, 20.8, 14.5, 8.3, and 3.9 
kJy for the PACSIOO, PACS160, SPIRE250, SPIRE350, 
and SPIRE500, respectively. For the LMC, these frac- 
tions are 0.91, 0.89, 0.87, 0.87, and 0.87 for global fluxes 
of 223, 259, 142, 73, and 31 kJy for the same bands. 
The global fluxes quoted here differ from those given by 
Meixner et al. (2013) due to our subtraction of MW cir- 
rus foreground and the additional background subtrac- 
tion step. 

The quantitative impact of correctly including the cor- 
related noise in the measurements can be illustrated by 
noting that assuming the noise is uncorrelated between 
bands results in the BEMBB model giving hts with a 
total SMC dust mass that is '^50% higher than the total 
dust mass given in Table 3. The importance of account- 
ing for the full likelihood function is equally important: 
the total SMC dust mass for the BEMBB model is ^50% 
higher using the ‘max’ values and ^30% lower using the 
‘exp’ values of log(E(i) when compared to the ‘realize’ 
value given in Table 3. The ‘realize’ values are the cor- 
rect values for determining the total dust mass values as 
they statistically reflect each pixel’s full likelihood func- 
tion, asymmetries and all, in the sum of the individual 
pixel masses. The ‘max’ and ‘exp’ values only reflect 
a limited portion of the likelihood function and this sys- 
tematically biases the results. This is additional evidence 
that the likelihood functions for are not well behaved 
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400 
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Figure 5. The gas-to-dust ratios (GDRs) are plotted as black 
circles for each of the three models and for both galaxies. The 
“reasonable” GDR range expected from scaling the MW diffuse to 
dense GDRs is given as a blue hatched region. The GDR range 
allowed by assuming the “maximum” depletions is given as a green 
hatched region (e.g. a lower limit on the GDR). 

Gaussians centered on the ‘max’ value (see Fig. 1). 

Our total dust masses are only lower limits as we do 
not include the dust responsible for the emission with 
surface brightnesses below 3a in any band. We can esti- 
mate the dust mass due to these <3cr regions by model- 
ing the integrated flux of these regions for each galaxy. 
Basically, we fit the SED that is the difference from the 
global fluxes quoted above and the integrated fluxes from 
<3cr pixels. The resulting integrated dust masses for the 
BEMBB model and the <3cr pixels are (5.9 ± 3.6) x 10^ 
and (1.6 ± 1.3) x 10^ Mq for the LMG and SMG, respec- 
tively. The uncertainties are quite large due to the low 
surface brightnesses and strong mixing of environments 
in these integrated SEDs. Gombining the <3a pixel dust 
masses with those for >3a pixel (Table 3), we find total 
dust masses of (7.3 ±1.7) x 10® and (8.3 ±2.1) x 10^ Mq 
for the LMG and SMG, respectively. For reference, the 
total gas masses that correspond to the same areas and 
same background removal as these total dust masses are 
3.1 X 10® and 3.0 x 10® Mq for the LMG and SMG, re- 
spectively. 

Bot et al. (2010a) obtained global dust masses for both 
galaxies by fitting Draine et al. (2007) dust models to 
their global fluxes. They found masses of 3.6 x 10® and 
0.29 — 1.1 X 10® Mq for the LMC and SMC, respectively. 
Leroy et al. (2007b) fit the spatially resolved Spitzer ob- 
servations with (Dale & Helou 2002) models and find a 
total SMC dust mass of 3 x 10® Mq. These values of 
the dust masses are factors of 4-5 larger than our values. 
The differences are likely due to different assumptions 
in the models used, the fitting techniques, the broader 
wavelength range of data, and/or the increased mixing 
of environments. 

6.3. Total Gas-to-Dust Ratios 

One test of the submm excess origin is to investigate 
how the overall gas-to-dust ratios for each model com- 
pare to the expected ratios. We explore overall gas-to- 
dust ratios as a test of the consistency of each dust model 
with expectations based on the measured gas masses and 
metallicities of the LMC and SMC. The detailed spatial 
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behavior of the gas-to-dust ratio with environment is in- 
vestigated in Roman-Duval et al. (this issue). 

The gas-to-dust ratios for each galaxy and all three 
models are given in Table 3. The dust masses are in- 
tegrated over all the pixels that are detected at >3cr 
in all observed bands. The total H gas masses given 
in the table footnote are integrated for the same pix- 
els as the dust masses. The HI masses are directly 
from the HI measurements (Stanimirovic et al. 2000; 
Muller et al. 2003; Kim et al. 2003) without any cor- 
rection for opaque HI (Dickey et al. 2000; Fukui et al. 
2014). The H2 masses are computed from CO ob- 
servations (Mizuno et al. 2001, 2006; Fukui et al. 2008; 
Wong et al. 2011) using Xco = 4.7 x 10^° (Hughes et al. 
2010) for the LMC and Xco = 6 x 10^^ (Leroy et al. 
2007a) for the SMC. The appropriate Xco to use is a 
matter of debate, but the expected range of this conver- 
sion factor is not large enough to strongly impact the to- 
tal gas masses (Fukui & Kawamura 2010; Bolatto et al. 
2013). The ratios given only include hydrogen, so are 
formally H gas-to-dust ratios, but for simplicity we refer 
to them as gas-to-dust ratios. 

The range of reasonable gas-to-dust ratios can be esti- 
mated three ways. The first scales the range of observed 
gas-to-dust ratios in the Milky Way by the LMC and 
SMC metallicities. The second assumes the Milky Way 
depletion factors and applies them to the measured LMC 
and SMC abundances. The third assumes all the met- 
als available are in the form of dust and this produces 
a minimum possible gas-to-dust ratio. The MW deple- 
tions and gas-to-dust ratios vary with environment and 
the global values in the Magellanic Clouds will be some 
unknown mix of different ISM environments. As a re- 
sult, we can only predict a possible range of gas-to-dust 
ratios. 

The first method assumes that the relative amount of 
metals in the LMC and SMC dust is the same as the 
MW, but scaled in proportion to each galaxy’s metal- 
licity. Thus, the expected gas-to-dust ratio will be 2X 
(LMC) and 5X (SMC) the MW gas-to-dust ratio. The 
MW gas-to-dust ratio varies from ~250 for the very dif- 
fuse ISM (F, = 0) to ~100 for the moderately dense 
ISM (F* = 1) (Jenkins 2009). For the LMC, we there- 
fore expect a gas-to-dust ratio between 200 to 500 while, 
for the SMC, we expect a gas-to-dust ratio between 500 
and 1250. The second method assumes the MW deple- 
tion patterns (Jenkins 2009) and the measured LMC and 
SMC abundances for each element (Russell & Dopita 
1992). The resulting expected LMC gas-to-dust ratios 
range between 150 to 360 and the expected SMC gas-to- 
dust ratios range between 540 to 1300. Combining the 
two different methods, the expected gas-to-dust ratios 
are 150 to 500 and 500 to 1300 for the LMC and SMC, 
respectively. Finally, the minimum allowed gas-to-dust 
ratio can be computed by assuming all the metals in the 
ISM in the form of dust. Assuming the measured LMC 
and SMC abundances, this gives minimum gas-to-dust 
ratios of 105 and 300, respectively. These expected gas- 
to-dust ratios are given in Table 3. 

The gas-to-dust ratios for all three models are plotted 
in Fig. 5 along with the allowed ranges for reasonable de- 
pletions and maximum depletion. From Table 3 and this 
figure, it is clear that the TTMBB models give gas-to- 


dust ratios that are lower than even possible assuming all 
the metals are present in dust. The TTMBB model gives 
low gas-to-dust ratios as it requires large dust masses for 
the second cold component to be able to reproduce the 
observed submm excess emission. Thus, the TTMBB 
model is not a reasonable model for the dust emission 
in the LMC or SMC. The SMBB and BEMBB models 
give similar gas-to-dust ratios for both galaxies. For the 
LMC, both models give ratios that are well within the 
reasonable range of values. For the SMC, these two mod- 
els both give values that are above the reasonable values. 
This is an indication that the depletions in the SMC are 
lower than those the in MW or that the dust proper- 
ties are different (e.g. a smaller Keff.ieo value than that 
assumed in this paper). 

6.4. Spatial Variations 

The spatial variations across both galaxies in the dif- 
ferent fit parameters for the BEMBB model are shown 
in Figs. 6 and 7. We only show the BEMBB results 
here as the evidence in the previous subsections gives a 
fairly strong indication that the BEMBB fits the data 
best (§6.1) and provides a physically reasonable gas-to- 
dust ratio (§6.3). The maps of dust surface density (E^) 
and temperature (Teff.d) show qualitatively similar be- 
haviors to previous works (Bot et al. 2004; Leroy et al. 
2007b; Bernard et al. 2008). In detail, our maps differ 
mainly in showing finer structure due to the higher spa- 
tial resolution Herschel observations. One illustration of 
this effect is that the peak Feff^d in the 30 Dor region in 
our map is ~60 K, significantly higher than the ^35 K 
found by Bernard et al. (2008). 

The higher spatial resolution of our maps does allow 
for detailed investigations of individual star forming re- 
gions. This is illustrated by Eig. 8 where cutouts of the 
BEMBB fit parameter maps for a star forming region in 
each galaxy are shown. The morphology of these two star 
forming regions is similar. The SPIRE 250 /im emission 
is strongly peaked in the region centers in contrast to the 
dust surface density which is more constant across the 
regions. This difference is caused by the center of these 
regions having high Teff^d values. The /3eff and 6500 maps 
of both regions have very similar morphologies, visually 
illustrating that these two fit parameters are strongly cor- 
related. Einally, the A{, images show coherent structures 
with fairly small variations overall. The submm excess as 
parametrized by 6500 is near zero in the center of the two 
star forming regions and rises rapidly to values around 
one near the edges. This behavior is intriguing, but the 
strong correlations of 6500 with /Jeff indicate that more 
work is needed to determine if this is real or due to noise 
induced correlations. 

The overall properties of the global submm excess be- 
tween the LMC and SMC show trends that are consistent 
with previous work. The average LMC and SMC 6500 
values are 0.27 and 0.43 when the average is done us- 
ing the ‘realize’ method and each pixel has equal weight. 
This can be visually seen in the 6500 images in Figs. 6 
and 7 where the SMC shows a higher filling factor of 
high 6500 values than the LMC. This trend of the lower 
metallicity SMC having a higher submm excess than the 
LMC is expected given the results from global studies 
of the submm excess Remy-Ruyer et al. (2013). A fairer 
comparison of the absolute value of 6500 with global SED 
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Figure 6. The spatial distribution of log(Eff), Teff^d? ^500» for BEMBB model are shown for the LMC using the 

expectation value for each pixel. In addition, the processed SPIRE 250 pm image (§2) is shown The images are shown using the cubehelix 
color mapping (Green 2011). The left/right and up/down streaks seen are residual instrumental artifacts that are aligned along the 
PACS/SPIRE scan direction. 
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Figure 7. The spatial distribution of Tefr,d» ^500? for the BEMBB model are shown for the SMC using the ‘exp’ value 

for each pixel. In addition, the processed SPIRE 250 fim image (§2) is shown The images are shown using the cubehelix color mapping 
(Green 2011). 
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Figure 8. The spatial distribution of Tefr,d» ^500? for the BEMBB model are shown for one star forming region each in 

the LMC and SMC using the ‘exp’ value for each pixel. In addition, the processed SPIRE 250 fim images (§2) is shown The images are 
shown using the cubehelix color mapping (Green 2011). 


fits is the dust surface density weighted averages that are 
0.11 and 0.26 for the LMC and SMC, respectively. Fi- 
nally, the average values of Af, are ~240 for both types 
of averages and both galaxies. This wavelength is simi- 
lar to that found by Li & Draine (2001) from fitting the 
DIRBE MW diffuse spectrum. 

To investigate the variations in fit parameters more 
quantitatively, we plot all the correlations between the 
different fit parameters for the LMC in Fig. 9. The plots 
for the SMC are very similar and are not shown. These 
plots show the density of points where each point rep- 
resents a single pixel. The values used for each pixel 
use the ‘realize’ method where the likelihood functions 
are randomly sampled once for each pixel. This means 
that these density plots statistically sample the full in- 
formation for the fit from each pixel. Repeating the ‘re- 
alize’ method process with a different random sampling 
for each pixel produces plots that are very similar. This 
indicates that these plots fully capture the correlations 
between fit parameters with a single sampling of each 
pixel’s likelihood function due to the large number of 
pixels. Plots created using the ‘max’ and ‘exp’ meth- 
ods are significantly different as they do not fully include 
the information on the uncertainties in the fits to each 
pixel. As an example of the difference between the differ- 
ent “best fit” methods, a flat likelihood function would 
show a single value for ‘max’ and ‘exp’, while the ‘real- 
ize’ method would have a value that was randomly dis- 
tributed over the entire parameter range. 

These plots show that many of the parameters are 
correlated with each other, sometimes quite strongly. 
The strongest correlations are seen between log(Sd) and 
Teff,d, log(Ed) and /leff.i, iPeff.d and /3eff,i, and ;5eff,i 
and 6500 . The origin of these correlations can be ei- 
ther real or a result of interactions between noise in the 
measurements and model fit parameters. The correla- 
tion between log(Sd) and Tes^d is real in that it reflects 
the detection thresholds of the HERITAGE data. Hot- 
ter dust can be detected a lower dust surface densities 


than cooler dust due to the ^ behavior of black- 
bodies. The anti-correlation between Teff,d and /3eff is 
one of the correlations that has been studied extensively 
to learn if it is due to noise or real variations in the 
dust properties (Dupac et al. 2003; Shetty et al. 2009a, b; 
Galliano et al. 2011; Juvela & Ysard 2012; Kelly et al. 
2012; Ysard et al. 2012; Veneziani et al. 2013). Labo- 
ratory data on dust analogs do show a shallow anti- 
correlation between Teff.d and PeS (Coupeaud et al. 
2011 ), but noise in measurements also produces a sim- 
ilar or larger anti-correlation. Kelly et al. (2012) have 
proposed to use a hierarchical Bayesian model to solve 
for the true TeS,d~PeS correlation, where the hierarchi- 
cal model assumes a single Tes^d and /3efi with some 
distribution around these values. In fitting an entire 
galaxy, such an assumption is not justified as, for exam- 
ple, there are regions near star formation that will be sig- 
nificantly hotter than regions further away. In addition, 
Juvela et al. (2013) hnd there are biases in all the cur- 
rently proposed methods for determining the true Teff.d” 
PeS relation. Thus, we choose to graphically display the 
correlations using the ‘realize’ method and not explicitly 
fit for the correlation. In future work, we plan to incor- 
porate additional observations of the ISM and physical 
models for the correlations between different ISM param- 
eters (e.g. dust and gas surface densities). 

Fig. 9 shows the correlations between the submm ex- 
cess 6500 and other dust properties. The value of 6500 
is positively correlated with and PeS.i and negatively 
correlated with Teff^d. This may be real or it may be due 
to the Teff^ versus PeS,i anti-correlation that is also very 
clearly seen. The positive correlation between 6500 and 
Ed is the opposite of what was found by Galliano et al. 
(2011) for a pathfinder study using a portion of the HER- 
ITAGE data on the LMG and Paradis et al. (2012) for 
the MW. The difference between these works and our 
work may be due to changes in the PACS and SPIRE cal- 
ibration, different fitting methods, and/or different dust 
emission models. Future work will investigate these dif- 
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Figure 9. The correlations for the LMC between all the five fit parameters for the BEMBB model are plotted. The plots are density 
plots where each point that contributes to the density is a single realization of the full likelihood function for a single pixel. 
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ferences by using the same data, same fitting code, and 
expanding the dust emission model to include more so- 
phisticated dust emission models. 

7. CONCLUSIONS AND FUTURE 

We find that the Magellanic Clouds show a submm ex- 
cess in the Herschel HERITAGE observations with a spa- 
tial resolution of ^10 pc. This submm excess seen in the 
Magellanic Clouds is more likely to be due to variations in 
the dust emissivity wavelength dependence than a second 
population of colder dust. This is based on the BEMBB 
model providing the best fit to the HERITAGE data and 
producing realistic gas-to-dust ratio values. The average 
submm excesses seen at 500 fim at ~10 pc resolution 
are 27% and 43% for the LMC and SMC, respectively. 
There are trends of the submm excess and environment 
(probed by and Tes^d), but the true nature of these 
trends will be investigated in future work incorporating 
more data and more physical models of the ISM. 

The total dust masses integrated over the pixels de- 
tected at 3cr in all five PACS/SPIRE bands using our 
favored model (BEMBB) are (7.3 ± 1.7) x 10® and 
(8.3 ± 2.1) X 10^ Mq for the LMC and SMC, respec- 
tively. These dust masses are significantly lower (fac- 
tors of 4-5) than would be expected from previous 
dust masses measurements (Leroy et al. 2007b; Bot et al. 
2010a). The lower dust masses we derive have impor- 
tant implications for the study of the lifecycle of dust in 
the Magellanic Clouds as the relative contributions be- 
tween Asymptotic Giant Branch (AGB), supernove, and 
the ISM for the formation of dust change significantly 
(Matsuura et al. 2009; Boyer et al. 2012; Matsuura et al. 
2013; Zhukovska & Henning 2013). 

Future work will focus on adding more physics to the 
fitting for dust properties. One rich area for future work 
will be to include constraints from other observations of 
the ISM in the Magellanic Clouds. An initial foray into 
this area is the focus of Roman-Duval et al. (this is- 
sue) who use the dust surface densities from this paper 
to investigate the dependence of the gas-to-dust ratio on 
environment. For the dust modeling in particular, future 
work will include more sophisticated dust grain models 
(e.g. Weingartner & Draine 2001; Compiegne et al. 2011; 
Galliano et al. 2011) and shorter wavelength infrared ob- 
servations (e.g., Spitzer IRAC/MIPS data) to better con- 
strain the possible grain compositions. 

We greatly benefited from conversations with Morgan 
Fouesneau, David Hogg, Derek Massa, and Daniel Weisz 
on the always interesting topic of fitting data with mod- 
els. We acknowledge financial support from the NASA 
Herschel Science Center, JPL contracts # 1381522 & 
1381650. M.R. acknowledges partial support from CON- 
ICYT project BASAL PFB-6. 

REFERENCES 

Aniano, G., Draine, B. T., Gordon, K. D., Sandstrom, K. 2011, 
PASP, 123, 1218 

Balog, Z., Muller, T., Nielbock, M., et al. 2013, Experimental 
Astronomy, arXiv:1309.6099 

Bendo, G. J., Griffin, M. J., Bock, J. J., et al. 2013, MNRAS, 

433, 3062 


Bernard, J.-P., Reach, W. T., Paradis, D., et al. 2008, AJ, 136, 
919 

Bianchi, S. 2013, A&A, 552, A89 

Boggess, N. W., Mather, J. C., Weiss, R., et al. 1992, ApJ, 397, 
420 

Bolatto, A. D., Wolfire, M., & Leroy, A. K. 2013, ARA&A, 51, 
207 

Bot, C., Boulanger, F., Lagache, G., Cambresy, L., & Egret, D. 
2004, A&A, 423, 567 

Bot, C., Ysard, N., Paradis, D., et al. 2010a, A&A, 523, A20 
Bot, C., Rubio, M., Boulanger, F., et al. 2010b, A&A, 524, A52 
Boudet, N., Mutschke, H., Nayral, G., et al. 2005, ApJ, 633, 272 
Boyer, M. L., Srinivasan, S., Riebel, D., et al. 2012, ApJ, 748, 40 
Cla37ton, G. C., Sz Martin, P. G. 1985, ApJ, 288, 558 
Compiegne, M., Verstraete, L., Jones, A., et al. 2011, A&A, 525, 
A103 

Coupeaud, A., Demyk, K., Meny, C., et al. 2011, A&A, 535, A124 
Dale, D. A., & Helou, G. 2002, ApJ, 576, 159 
Dale, D. A., Aniano, G., Engelbracht, C. W., et al. 2012, ApJ, 
745, 95 

Desert, F.-X., Boulanger, F., & Puget, J. L. 1990, A&A, 237, 215 
Dickey, J. M., Mebold, U., Stanimirovic, S., & Staveley-Smith, L. 
2000, ApJ, 536, 756 

Draine, B. T., Li, A. 2007, ApJ, 657, 810 

Draine, B. T., Dale, D. A., Bendo, G., et al. 2007, ApJ, 663, 866 

Draine, B. T., Aniano, G., Krause, O., et al. 2014, ApJ, 780, 172 

Dupac, X., Bernard, J.-P., Boudet, N., et al. 2003, A&A, 404, Lll 

Fitzpatrick, E. L. 1985, ApJ, 299, 219 

Fukui, Y., &; Kawamura, A. 2010, ARA&A, 48, 547 

Fukui, Y., Kawamura, A., Minamidani, T., et al. 2008, ApJS, 

178, 56 

Fukui, Y., Okamoto, R., Yamamoto, H., et al. 2014, ArXiv 
e-prints, arXiv: 1401.7398 

Galametz, M., Madden, S. C., Galliano, F., et al. 2011, A&A, 

532, A56 

Galametz, M., Hony, S., Galliano, F., et al. 2013, MNRAS, 431, 
1596 

Galliano, F., Madden, S. C., Jones, A. P., Wilson, C. D., Sz 
Bernard, J.-P. 2005, A&A, 434, 867 
Galliano, F., Madden, S. C., Jones, A. P., et al. 2003, A&A, 407, 
159 

Galliano, F., Hony, S., Bernard, J.-P., et al. 2011, A&A, 536, A88 
Gordon, K. D., & Clayton, G. C. 1998, ApJ, 500, 816 
Gordon, K. D., Clayton, G. C., Misselt, K. A., Landolt, A. U., & 
Wolff, M. J. 2003, ApJ, 594, 279 

Gordon, K. D., Galliano, F., Hony, S., et al. 2010, A&A, 518, L89 
Gordon, K. D., Meixner, M., Meade, M. R., et al. 2011, AJ, 142, 
102 

Green, D. A. 2011, Bulletin of the Astronomical Society of India, 
39, 289 

Griffin, M. J., Abergel, A., Abreu, A., et al. 2010, A&A, 518, L3 
Griffin, M. J., North, C. E., Schulz, B., et al. 2013, MNRAS, 434, 
992 

Gut, A. 2009, An Intermediate Course in Probability, Springer 
Texts in Statistics (Springer) 

Herschel Space Observatory. 2011, SPIRE Observers Manual, 
Tech. Rep. HERSCHEL-DOC-0798, Version 2.4, ESA, 
Noordwijk 

Hildebrand, R. H. 1983, QJRAS, 24, 267 

Hilditch, R. W., Howarth, I. D., & Harries, T. J. 2005, MNRAS, 
357, 304 

Hughes, A., Wong, T., Ott, J., et al. 2010, MNRAS, 406, 2065 
Israel, F. P., Wall, W. F., Raban, D., et al. 2010, A&A, 519, A67 
Jager, C., Mutschke, H., & Henning, T. 1998, A&A, 332, 291 
Jenkins, E. B. 2009, ApJ, 700, 1299 

Juvela, M., Montillaud, J., Ysard, N., Sz Lunttila, T. 2013, A&A, 
556, A63 

Juvela, M., Ysard, N. 2012, A&A, 541, A33 
Kelly, B. C., Shetty, R., Stutz, A. M., et al. 2012, ApJ, 752, 55 
Kessler, M. F., Steinz, J. A., Anderegg, M. E., et al. 1996, A&A, 
315, L27 

Kim, S., Staveley-Smith, L., Dopita, M. A., et al. 2003, ApJS, 

148, 473 

Kirkpatrick, A., Calzetti, D., Galametz, M., et al. 2013, ApJ, 778, 
51 

Lequeux, J., Maurice, E., Prevot-Burnichon, M.-L., Prevot, L., &; 
Rocca-Volmerange, B. 1982, A&A, 113, L15 


20 


Gordon et al. 


Leroy, A., Bolatto, A., Stanimirovic, S., et al. 2007a, ApJ, 658, 
1027 

Leroy, A., Cannon, J., Walter, F., Bolatto, A., &: Weiss, A. 

2007b, ApJ, 663, 990 

Li, A., & Draine, B. T. 2001, ApJ, 554, 778 
Maiz Apellaniz, J., & Rubio, M. 2012, A&A, 541, A54 
Mather, J. C., Fixsen, D. J., & Shafer, R. A. 1993, in Society of 
Photo-Optical Instrumentation Engineers (SPIE) Conference 
Series, Vol. 2019, Infrared Spaceborne Remote Sensing, ed. 

M. S. Scholl, 168M79 

Mathis, J. S., Rumpl, W., & Nordsieck, K. H. 1977, ApJ, 217, 425 
Matsuura, M., Woods, P. M., & Owen, P. J. 2013, MNRAS, 429, 
2527 

Matsuura, M., Barlow, M. J., Zijlstra, A. A., et al. 2009, 

MNRAS, 396, 918 

Meixner, M., Gordon, K. D., Indebetouw, R., et al. 2006, AJ, 132, 
2268 

Meixner, M., Galliano, F., Hony, S., et al. 2010, A&A, 518, L71 
Meixner, M., Panuzzo, P., Roman-Duval, J., et al. 2013, AJ, 146, 
62 

Mennella, V., Brucato, J. R., Colangeli, L., et al. 1998, ApJ, 496, 
1058 

Mennella, V., Colangeli, L., Bussoletti, E. 1995, A&A, 295, 165 
Misselt, K. A., Clayton, G. G., & Gordon, K. D. 1999, ApJ, 515, 
128 

Mizuno, N., Muller, E., Maeda, H., et al. 2006, ApJ, 643, L107 
Mizuno, N., Rubio, M., Mizuno, A., et al. 2001, PASJ, 53, L45 
Muller, E., Staveley-Smith, L., Zealey, W., Sz Stanimirovic, S. 
2003, MNRAS, 339, 105 

Muller, T., Nielbock, M., Balog, Z., Klaas, U., Sz Vilenius, E. 
2011a, PAGS Photometer - Point-Source Flux Calibration, 

Tech. Rep. PICC-ME-TN-037, Herschel 
Muller, T., Okumura, K., Sz Klaas, U. 2011b, PACS Photometer 
Passbands and Colour Correction Factors for Various Source 
SEDs, Tech. Rep. PICC-ME-TN-038, Herschel 
Paradis, D., Bernard, J., Sz Meny, C. 2009, A&A, 506, 745 
Paradis, D., Paladini, R., Noriega-Crespo, A., et al. 2012, A&A, 
537, A113 

Peimbert, A., Sz Peimbert, M. 2010, ApJ, 724, 791 
Pilbratt, G. L., Riedinger, J. R., Passvogel, T., et al. 2010, A&A, 
518, LI 

Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2011, 
A&A, 536, A17 

Poglitsch, A., Waelkens, C., Geis, N., et al. 2010, A&:A, 518, L2 


Prevot, M. L., Lequeux, J., Prevot, L., Maurice, E., Sz 
Rocca-Volmerange, B. 1984, ASzA, 132, 389 
Reach, W. T., Dwek, E., Fixsen, D. J., et al. 1995, ApJ, 451, 188 
Remy-Ruyer, A., Madden, S. G., Galliano, F., et al. 2013, A&A, 
557, A95 

Roth, K. G., Sz Blades, J. C. 1997, ApJ, 474, L95 
Russell, S. C., Sz Dopita, M. A. 1992, ApJ, 384, 508 
Schwering, P. B. W. 1989, ASzAS, 79, 105 
Schwering, P. B. W., Sz Israel, F. P. 1989, A&;AS, 79, 79 
Shetty, R., Kauffmann, J., Schnee, S., Sz Goodman, A. A. 2009a, 
ApJ, 696, 676 

Shetty, R., Kauffmann, J., Schnee, S., Goodman, A. A., Sz 
Ercolano, B. 2009b, ApJ, 696, 2234 
Skibba, R. A., Engelbracht, C. W., Aniano, G., et al. 2012, ApJ, 
761, 42 

Sofia, U. J., Gordon, K. D., Cla 3 don, G. C., et al. 2006, ApJ, 636, 
753 

Stanimirovic, S., Staveley-Smith, L., van der Hulst, J. M., et al. 
2000, MNRAS, 315, 791 

Staveley-Smith, L., Kim, S., Calabretta, M. R., Haynes, R. F., Si 
Kesteven, M. J. 2003, MNRAS, 339, 87 
Valencic, L. A., Cla 3 don, G. C., Sz Gordon, K. D. 2004, ApJ, 616, 
912 

Valencic, L. A., Clayton, G. C., Gordon, K. D., Sz Smith, T. L. 
2003, ApJ, 598, 369 

Veneziani, M., Piacentini, F., Noriega-Crespo, A., et al. 2013, 
ApJ, 772, 56 

Walker, A. R. 2012, Ap&SS, 341, 43 
Weingartner, J. C., Sz Draine, B. T. 2001, ApJ, 548, 296 
Welty, D. E., Sz Crowther, P. A. 2010, MNRAS, 404, 1321 
Welty, D. E., Lauroesch, J. T., Blades, J. C., Hobbs, L. M., Sz 
York, D. G. 1997, ApJ, 489, 672 
— . 2001, ApJ, 554, L75 

Werner, M. W., Roellig, T. L., Low, F. J., et al. 2004, ApJS, 154, 
1 

Wong, T., Hughes, A., Ott, J., et al. 2011, ApJS, 197, 16 
Wright, E. L., Mather, J. C., Bennett, C. L., et al. 1991, ApJ, 
381, 200 

Ysard, N., Juvela, M., Demyk, K., et al. 2012, ASzA, 542, A21 
Zhukovska, S., Sz Henning, T. 2013, A&:A, 555, A99 
Zubko, V., Dwek, E., &; Arendt, R. G. 2004, ApJS, 152, 211 
Zubko, V. G., Mennella, V., Colangeli, L., & Bussoletti, E. 1996, 
MNRAS, 282, 1321 


